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+-• Abstract 

P i 

Q We study the formation of inhomogeneous order in the Einstein-Maxweh- 

,-^ axion system, dual to a 2+1 dimensional field theory that exhibits a sponta- 

neously generated current density, momentum density and modulated scalar 
operator. Below the critical temperature, the Reissner-Nordstrom-AdS black 
hole becomes unstable and stripes form in the bulk and on the boundary. The 

(T^ bulk geometry possesses striking geometrical features, including a modulated 

horizon that tends to pinch off as T — )■ 0. On a domain of fixed length, we find 
a second order phase transition to the striped solution in each of the grand 
canonical, canonical and microcanonical ensembles, with modulated charges 
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(y-^ that grow and saturate as we lower the temperature and descend into the in- 

T-H homogeneous phase. For the black hole on an infinite domain, a similar second 

^ order transition occurs, and the width of the dominant stripe increases in the 

*k> zero temperature limit. 



1 Introduction and Summary 

The gauge-gravity duality is a relationship between a strongly coupled field theory 
and a gravity system in one higher dimension. This correspondence has been fruitful 
in studying various field theory phenomena by translating the problem to the gravi- 
tational context. In particular, the duality has shone new light on many condensed 
matter systems - see [H El |3l H] for reviews. 

Early models in this area, such as the holographic superconductor [5], focused 
on homogeneous phases of field theories. In this case, the fields on the gravity side 
depend only on the radial coordinate in the bulk and the problem reduces to the 
solution of ODEs. However, many interesting phenomena occur in less symmetric 
situations. Generically, the problem of finding the gravity dual to an inhomogeneous 
boundary system will necessitate solving relatively more difficult PDEs, almost al- 
ways resulting in the need for numerical methods. While these become technically 
hard problems, there exist established numerical approaches. Due to the success 
of the holographic method in studying homogeneous situations, it is worthwhile to 
push the correspondence to these less symmetric situations in order to describe more 
general phenomena in this context. 

One particular area of condensed matter that appears to be amenable to a holo- 
graphic description is the appearance of striped phases in certain materials]^ These 
phases are characterized by the spontaneous breaking of translational invariance 
in the system. Examples include charge density waves and spin density waves in 
strongly correlated electron systems, where either the charge and/or the spin densi- 
ties become spatially modulated (for a review see [8]). The formation of stripes is 
conjectured to be related to the mechanism of superconductivity in the cuprates ^. 
To approach this striking phenomenon from the holographic perspective, one would 
look for an asymptotically AdS gravity system which allows a spontaneous transition 
to a modulated phase. 

Recently, several interesting spatially modulated holographic systems have been 
studied. One way to study stripes on the boundary is to source them by imposing 
spatial modulation in the non-normalizable modes of some fields, explicitly breaking 
the translation invariance, as in (TUl [II]n However, if one wishes to make contact 
with the context described above, it is important that the inhomogeneity emerges 
spontaneously rather than be introduced explicitly. 

In some cases, the spatially modulated phase has an extra symmetry, allowing 



^Stripes are also known to form in large N QCD [6, 7 . 

^In a similar vein, more recently, lattice-deformed black branes have been of interest in studies 
of conductivity in holographic models [THl ESI Ell IH] • 



the situation to be posed as a co-homogeneity one problem on the gravity side. 
Examples include systems in which one of the translational Killing vectors is replaced 
by a helical Killing vector [161 [13 UHl [19], EOl El]- More general inhomogeneous 
instabilities, in which one of the translation symmetries is fully broken, have been 
described in a phenomenological model [22] and in certain H^ND = 6 brane systems 



In this work, we study the full non-linear co-homogeneity two striped solutions to 
the Einstein-Maxwell-axion model that stem from the normalizable, inhomogeneous 
modes of the Reissner-Nordstrom-AdS solution detailed in [22] . In this model, below 
a critical temperature, stripes spontaneously form in the bulk and on the boundary. 
We study the properties of the stripes in both the fixed length system, in which the 
wavenumber is set by the size of the domain and charges are integrated over the 
stripe, and the infinite system, in which the corresponding thermodynamic densities 
are studied. For the black hole at fixed length, we examine the behavior in different 
thermodynamic ensembles as we vary the temperature and wavenumber. 

The study is facilitated by a numerical solution to the set of coupled Einstein 
and matter equations in the bulk. Inspired by the black string case [SUl EI], we fix 
the metric in the conformal gauge, resulting in a set of field equations and a set of 
constraint equations. Then, as described in [30], the resulting constraint equations 
can be solved by imposing particular boundary conditions on the fields. 

As well as being of interest from the holographic perspective these numerical 
solutions are important as they represent new inhomogeneous black hole solutions 
in Ads. We find strong evidence that the unstable homogeneous branes transition 
smoothly to the striped state below the critical temperature]^ As we approach zero 
temperature the relative inhomogeneity is seen to grow without bound and the black 
hole horizon tends to pinch off, signalling the formation of a spacetime singularity in 
this limit. 

A subset of our results has already been announced in [33], in this paper we 
provide full details. The summary of the results follow: 

Boundary field theory 

• We calculate the fully back-reacted normalizable inhomogeneous modes. 

• The stripes have momentum, electric current and modulations in charge and 
mass density (see [36] for a recent study of angular momentum generation). 



•^ Other studies of inhomogeneity in the context of holography include jSHl [13 [251l29j . 

* The instability to the formation of the striped black branes resembles the black string instability 



which is known to be of the second order for high enough dimensions 



• As a function of temperature, the modulations start small, then grow and 
saturate as T — )■ 0. 

• We study the stripe of fixed length in various ensembles, finding a second or- 
der phase transition, for sufficiently large axion coupling, in each of the grand 
canonical (temperature T, chemical potential /j, fixed), canonical (T, charge N 
fixed) and microcanonical (mass M, A^ fixed) ensembles. We compute corre- 
sponding critical exponents. 

• For the infinite length system, there is a second order transition to a striped 
phase. The width of the dominant stripe grows as the temperature is decreased. 

• In the zero temperature limit, within the accuracy of our numerics, the entropy 
appears to approach a non-zero value. 

Bulk geometry 

The new inhomogeneous black brane solutions that we find have peculiar features, 
including 

• The inhomogeneities are localized near the horizon, and die off asymptotically 
following a power law decay. 

• The phenomena of vorticity, frame dragging and the magneto-electric effect 
similar to one produced by a near horizon topological insulator are observed. 

• The inhomogeneous black brane has a neck and a bulge. In the curvature at 
the horizon, the maximum is at the bulge. In the limit of small temperatures, 
the neck shrinks to zero size. 

• The proper length of the horizon grows when temperature is decreasing, and 
diverges as l/T"'^ in the limit T — ?► 0. The proper length in the stripe direction 
increases from the boundary to the horizon, which can be thought of as a 
manifestation of an "Archimedes effect" . 

In [|2| we define our model and set up our numerical approach, describing our 
ansatz, boundary conditions and solving procedure. Then, in ^ we report on in- 
teresting geometrical features of the bulk solutions, ^studies the solutions at fixed 
length from the point of view of the boundary theory. There, we make the compar- 
ison to the homogeneous solution and find a second order transition, in addition to 
describing the observables in the theory. In ^ we relax the fixed length condition 



and find the striped solution that dominates the thermodynamics for the infinite 
system. Appendix |A] provides details about computing the observables of the inho- 
mogeneous solutions while appendix |B] gives more details on the numerics, including 
checks of the solutions and validations of our numerical method. 

Note added: As this manuscript was being completed, [37] and [38] appeared, 
which use a different method and have some overlap with this work. 

2 Numerical set-up: Einstein- Maxwell-axion model 

In [22], perturbative instabilities of the Reissner-Nordstrom-AdS (RN for short) black 
brane were found within the Einstein-Maxwell-axion model. In [33] and here, we 
construct the full non-linear branch of stationary solutions following this zero mode. 

2.1 The model and ansatz 

The Lagrangian describing our coupled system can be written as [22] 

C = 1{R + 12) - ^d'^^pd.^P - ^mV^ - \f^'F,^ - J^ ^^'^ i, 6^-^-F,.F,., (2.1) 

where R is the Ricci scalar, F^^ is the Faraday tensor, ip is a pseudo-scalar field and 
g is the determinant of the metric. We use units in which the AdS radius P = 1/2, 
Newton's constant 8ttGn = 1 and c = h = 1, and choose m^ = —4. The constant ci 
controls the strength of the axion coupling. 

For this choice of scalar field mass, instabilities exist for all choices of Ci. For 
Ci = 0, the instability is towards a black hole with neutral scalar hair. For ci > 
0, inhomogeneous instabilities along one field theory direction exist for a range of 
wavenumbers k. The critical temperature at which each mode becomes unstable 
depends on the wavenumber: Tc{k). For a given Ci, there is a maximum critical 
temperature, above which there are no unstable modes. As one increases ci, the 
critical temperature of a given mode k increases, such that for a fixed temperature 



a larger range of wavenumbers will be unstable. See appendix B.l for more details 
on the perturbative analysis. 

One may consider generalizations of this action, including higher order couplings 
between the scalar field and the gauge field. In particular, as discussed in p2j . 
generalizing the Maxwell term as —^-^F^'^F^,^, where T{ip) is a function of the scalar 
field, results in a model that can be uplifted to a D = 11 supergravity solution (for 
particular choices of ci, m, and the parameters in T{tp)). In this study, we wish to 



study the formation of holographic stripes phenomenologically. The existence of the 
axion-couphng term (ci 7^ 0) is a sufficient condition for the inhomogeneous solutions 
and so we set T{tp) = 1 here. 

We are looking for stationary black hole solutions that can be described by an 
ansatz of the form 

ip = ilj{r,x), A = At{r,x)dt + Ay{r,x)dy, (2.2) 

where r is the radial direction in AdS and x is the field theory direction along which 
inhomogeneities form. We term the scalar field and gauge fields collectively as the 
matter fields, /(r) is a given function whose zero defines the black brane horizon. 
We take /(r) to be that of the RN solution, 

so that the horizon is located at r = tq. The homogeneous solution is the RN black 
brane, given by 

A = B = C = W = ij = Ay = 0, A{r)=fi{l-ro/r), (2.4) 

where /i is the chemical potential. Above the maximum critical temperature, this is 
the only solution to the system. 

To find the non-linear inhomogeneous solutions, we numerically solve the equa- 



tions of motion derived from the ansatz (2.2). The Einstein equation results in four 
second order elliptic equations, formed from combinations of Gj— T/ = 0, Gy—Ty = 0, 
G^^ — Ty = 0, and G^ + G^ — (T^ + T^) = 0, and two hyperbolic constraint equations, 
G^ — T^ = and G^. — G^ — {T^' — T^) = 0, for the metric functions. The gauge 
field equations and scalar field equation give second order elliptic equations for the 



matter fields. For completeness, the full equations are given in appendix B.2 Our 
strategy will be to solve these seven elliptic equations subject to boundary conditions 
that ensure that the constraint equations will be satisfied on a solution. Below, we 
describe the constraint system and our boundary conditions. For more details about 
the numerical approach, we refer to appendix |B} 

2.2 The constraints 

The two equations G^ — T^ = and G^. — G^ — (T^ ~ T^) = 0, which we do not 
explicitly solve, are the constraint equations. Using the Bianchi identities [3D], we 



see that the constraints satisfy 

d. {v^{G: - TD) + 2r2 v/79, (r^ v/7^/=^(G; - G^ - (X; - T;))) = 0, (2.5) 



2r^ V7a. (v/=^(G^. - t;)) - a. (r^ v7v^(g; - g^ - (t; - t;))) = 0. (2.6) 

Defining f hy df = 2r^^dr gives Cauchy-Riemann relations 

d. {V^iGl - T:)) + dr [r'^ffV^giC: -Gl~ (T; - T^))) = 0, (2.7) 



5. (v^lG^ - T:)) - d, (r'^rfV^g{G: -Gl- (T; - T^))) = 0, (2.8) 

showing that the weighted constraints satisfy Laplace equations. Then, satisfying 
one constraint on the entire boundary and the other at one point on the boundary 
implies that they will both vanish on the entire domain. In practice we will take either 
zero data or Neumann boundary conditions at the boundaries in the x-direction. The 
unique solution to Laplace's equation with zero data on the horizon and the boundary 
at infinity and these conditions in the x-direction is zero. Therefore, as long as we 
fulfill one constraint at the horizon and the asymptotic boundary and the other at 
one point (on the horizon or boundary), the constraints will be satisfied if the elliptic 
equations are. Our boundary conditions will be such that a/^(G^ — T^) = at the 
horizon and conformal infinity and that r'^y/Jy/—g{G'r — G^ — {T^ — T^)) = at one 
point on the horizon. 

2.3 Boundary conditions 

The elliptic equations to be solved are subject to physical boundary conditions. 
There are four boundaries of our domain (see Fig. fl]): the horizon, the conformal 
boundary, and the periodic boundaries in the x-direction, which are described next. 

Staggered periodicity 

To specify the boundary conditions in the x direction we look at the form of the 



linearized perturbation which becomes unstable (see appendix B.l ). To leading order 
in the perturbation parameter A, they are of the form: 

ip{x) ~ Acos(/cx), 
Ay{x) ~ Asin(A;x), 
gty{x) ~Asin(A;x), (2.9) 
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Figure 1: A summary of the boundary conditions on our domain. At the horizon, 
r = ro, we impose regularity conditions. At the conformal boundary, r — )■ oo, we 
have fall off conditions on the fields (imposed at large but finite r = Tcut) such that we 
do not source the inhomogeneity. In the x-direction, we use symmetries to reduce the 
domain to a quarter period L/4. Then, we impose either periodic or zero conditions 
on the fields, according to their behavior under the discrete symmetries discussed in 
the text, {h collectively denotes the fields {gtt,gxx,gyy,At}.) In addition to these, 
we explicitly satisfy the constraint equation y/—gGl. = on the horizon and the 
conformal boundary. 



where k is the wavenumber of the unstable mode. To second order in the perturbation 
parameter, the functions gu,gxx,gyy and At (which we denote collectively as h) are 
turned on, with the schematic behavior 



h{x) r^ X^{cos{2kx) + C), 



(2.10) 



where C are independent of x. 

All these functions are periodic with period L = ^. However, they are not 
the most general periodic functions with period L. For numerical stability it is 
worthwhile to specify their properties further and encode those properties in the 
boundary conditions we impose on the full solution. We concentrate on the behavior 
of the perturbation with respect to two independent Z2 reflection symmetries. 



The first Z^ symmetry is that of x — )■ — x, y -^ — y, which is a rotation in the x, y 
plane. This is a symmetry of the action and of the hnearized perturbation (keeping in 
mind that Ay and gty change sign under reflection of the y coordinate). We conclude 
therefore that this is a symmetry of the full solution. 

Similarly, the Z2 operation x — )■ ^ — x, y — )■ —y is a symmetry of the action, 
which is also a symmetry of the linearized system when accompanied by A — )■ —A. 
In other words the functions V', Ay, gty are restricted to be odd with respect to this 
Z2 operation, while the rest of the functions, which we collectively denoted as /i, are 
even. 

The two symmetries defined here restrict the form of the functions that can 
appear in the perturbative expansions for each of the functions above. For example, 
it is easy to see that the function '^{x) gets corrected only in odd powers of A and the 
most general form of the harmonic that can appear in the perturbative expansion is 
cos(?7./cx), for n odd. Similar comments apply to the other functions above. 

We restrict ourselves to those harmonics which may appear in the full solution. 
The most efficient way to do so is to work with a quarter of the full period L 
(reconstructing the full periodic solution using the known behavior of each function 
with respect to the two Z^ operations defined above). The specific properties of 
each function appearing in our solutions are imposed by demanding the following 
boundary conditions: 





d^i\}[x = 0) 


= 0, 


Hx = ^)=0 




Ay{X = 0) 


= 0, 


d^Ay{x = -) 




gty{x = 0) 


= 0, 


dxgtyix = -) 




d^x = 0) 


= 0, 


dMx = 4 ) = 


At the horizon 









(2.11) 



In our coordinates (2.2) the horizon is at fixed r = tq. For numerical convenience 
we introduce another radial coordinate p = -^/r^ — Tq, such that the horizon is at 
p = Ojj Expanding the equations of motion around p = yields a set of Neumann 
regularity conditions, 

d,A = d,C = dpW = d,ilj = dpAt = dpAy = 0, (2.12) 



^In the rest of the paper, we use r and p interchangeably as our radial coordinate. We use the 
coordinate p in the numerics. 

9 



and two conditions in the inhomogeneous direction along the horizon, 

d,W = d,{At + WAy) = {}. (2.13) 

Thus, both W and the combination At + W Ay are constant along the horizon. The 



boundary conditions in the x direction (2.11) imply that ly = 0. Then, the second 



condition together with regularity of the vector field A on the Euclidean section give 
that ^t = on the horizon. 

The regularity conditions give eight conditions for the six functions A, C,W,ijj, At 
and Ay. In principle, we would choose any six of these to impose at the horizon. If we 
find a non-singular solution to the equations, then the other two conditions should 
also be satisfied. In practice, some of these conditions work better than others for 
finding the numerical solution. We find that using Neumann conditions for A, C, ip, 
and A„ and Dirichlet conditions for W and At results in a more stable relaxation'^ 



The conditions for B are determined using the constraint equations. Expanding 
the weighted constraints at the horizon, we find 

v^(G;-T;) (X d^{A-B) + Oip), (2.14) 

r'v7v^(G::-G^-(r;-T;)) ex d,B + o{p). (2.15) 

The first condition gives constant surface gravity (or temperature) along the horizon. 
As discussed above, we will impose one constraint at the horizon and the boundary, 
and the other at one point. In practice, we will satisfy r"^ \fj y/^iG^^—G^ — {T^ —T^)) 
at (p, x) = (0, 0), updating the value of B at this point using the Neumann condition 
dpB = 0. This will set the difference {B — A)|(p,x)=(o,o) = do, which we will then 
use to update B using a Dirichlet condition along the rest of the horizon, satisfying 

V^iGi - r;) = 0. 

At the conformal boundary 

In our coordinates, the boundary is at r = oo. Since we are looking for spontaneous 
breaking of homogeneities, our boundary conditions will be such that the field theory 
sources are homogeneous. This implies that the non-normalizable modes of the bulk 
fields are homogeneous. The inhomogeneity of the striped solutions will be imprinted 
on the normalizable modes of the fields, or the coefficient of the next-to-leading fall- 
off term in the asymptotic expansions. 



^Using Neumann conditions at the horizon for W and At resuhs in values at the horizon that 
converge to zero with step-size, consistent with the above analysis. 



10 



The form of our metric ansatz is such that the metric functions A,B,C and 
W represent the normahzable modes of the metric. Imposing that the geometry is 
asymptotically AdS with Minkowski space on the boundary implies that these four 
metric perturbations must vanish as r — )■ oo. By expanding the equations of motion 
near the boundary, one can show that A,B,C and W fall off as 1/r^. In practice, 
we place the outer boundary of our domain at large but finite rcut and impose the 
fall-off conditions there. 

As in the RN solution, we source the field theory charge density with a homoge- 
neous chemical potential, corresponding to a Dirichlet condition for the gauge field 
At at the boundary. In the inhomogeneous solutions, we expect the spontaneous 
generation of a modulated field theory current jy{x), dual to the normahzable mode 
of Ay. Solving the equations near the boundary with these conditions reveals the 
expansions At = fi + 0{l/r) and Ay = 0(l/r), which we impose numerically at Vcut- 

The scalar field equation of motion gives the asymptotic solution 

^ = 7xr + ^ + ---' (2-16) 



where 

A 



:^(3±v/9 + 4(/m)2). (2.17) 

For the range of scalar field masses —9/2 < m^ < —5/2, both modes are normahz- 
able, and fixing one mode gives a source for the other. In our study we will choose 
m^ = —4, giving A_ = 1, A+ = 2. Since we are looking for spontaneous symmetry 
breaking, in this case we must choose either ^^^^ = or ■ip^'^^ = 0. We choose the 
former, so that ip falls off as 1/r^. 

Now, consider the weighted constraint ^—gG''^. As discussed above, in order to 
solve the constraint system, we require this to disappear at the conformal boundary. 
Near the boundary, a/— ^ oc r^ + . . . , so for y/—gGl to disappear we must have 
Gl = 0(l/r^). Expanding the equations near the boundary we have 

where X = X^^\x)/r^ + ... for X = {A, i?,C}. Therefore, in addition to the 
boundary conditions mentioned above, for ^f—gG"^^ = to be satisfied at r = 00, 
it appears that we should have that 'iA^'^\x) + 2B^^\x) + 3C*-^''(x) = const. The 



means to impose this addition condition comes from the fact that our metric (2.2) 
has an unfixed residual gauge freedom [SH], allowing one to transform to new f = 
f{r, x),x = x{r, x) coordinates which are harmonic functions of r and x. Performing 

11 



such a transformation generates an additional function in (2.18), which can then be 
chosen to ensure that the constraint is satisfied (in appendix |B| we describe how). 
This condition imphes the conservation of the boundary energy momentum tensor, 
see appendix \K\ 

2.4 Parameters and algorithm 

The physical data specifying each solution is the chemical potential /x, the temper- 
ature T, and the periodicity Lj^ Since the boundary theory is conformal, it will 
only depend on dimensionless ratios of these parameters. This manifests itself in the 
following scaling symmetry of the equations: 

r^Xr, {t,x,y) ^ j{t,x,y), A^ ^ XA,,. (2.19) 

We use this to select /z = 1. Then, our results are functions of the dimensionless 
temperature T/fi and the dimensionless periodicity Lfi. 

The temperature is controlled by the coordinate location of the horizon. For a 
given To, the temperature of the RN phase is Tq = (l/87rro)(12rQ — 1) while the tem- 
perature of the inhomogeneous solution is T = e~°'°To. Recall that {B — A) \rg = rfg is 
dynamically generated by satisfying the constraints at the horizon. From our numer- 
ical solutions, we find that do monotonically increases as we lower the temperature, 
so that To gives a reliable parametrization of the physical temperature T. In practice, 
we generate solutions by choosing values of Tq below the critical temperature Tc{k). 

We solve the equations by finite-difference approximation (FDA) techniques. We 



use second order FDA on the equations (B.6) - (B.12) before using a point-wise 
Gauss-Seidel relaxation method on the resulting algebraic equations. For the results 
in this paper, for ci = 4.5, a cutoff of pcut = {6, 8} was used while for ci = 5.5 and 
Ci = 8, for which the modulations were larger, a cutoffs of pcut = 10 and pcut = 12 
correspondingly were used. Grid spacings used for the FDA scheme were in the range 
dp, dx = 0.04 — 0.005. Neumann boundary conditions are differenced to second order 
using one-sided FDA stencils in order to update the boundary values at each step. 
At the asymptotic boundary pcut we impose the boundary conditions by second order 
differencing a differential equation based on the fall-off (for example, drA = —3A/r) 
to obtain an update rule for the boundary value. As a result we find quadratic 
convergence as a function of grid-spacing for our method, see appendix |B. 5.1 



^Fixing /i, T and L gives the system in the grand canonical ensemble. Once the phase space has 
been mapped in one ensemble other ensembles can be considered via appropriate reinterpretation 
of the numerical data. See 34| for a description of this process. 
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3 The solutions 



The system of equations (B.6 B.12 ) is solved subject to boundary condition described 
in the previous sections. The details of our numerical algorithm are found in appendix 
[B| Here we focus on the properties of the solutions and their geometry. 

Unless otherwise specified the following plots were obtained using the axion cou- 
pling of c\ = 4.5. In this section, we consider solutions for which the periodicity 
is determined by the dominant critical wavenumber k^, for ci = 4.5, this gives 
L/i/4 ~ 2.08, see Table [I] We found that the geometry and most of the other 
features are qualitatively similar for the couplings Ci = 5.5 and Ci = 8. A convenient 
way to parametrize our inhomogeneous solutions is by the dimensionless temperature 
T/Tc, relative to the critical temperature Tc, below which the translation invariance 
along X is broken. For ci = 4.5, our method allows us find solutions in the range 
0.003 < T/T, < 0.9. 

3.1 Metric and fields 

For subcritical temperatures, as we descend into inhomogeneous regime, the metric 
and the matter fields start developing increasing variation in x. Fig. |2] displays 
the metric functions, and Fig. [3] shows the non vanishing components of the vector 
potential field and of the scalar field for T/Tc ^ 0.11 over a full period in the x 
direction. The variation of all fields is maximal near the horizon of the black hole 
at p = Y r^ — Tq = 0, and it gradually decreases toward the conformal boundary, 
p — > oo. 

Many of the special features of the solutions we find may be explained via axion 
electrodynamics as seen in the effective description of the electromagnetic response 
of a topological insulator. This effect is mediated by the interaction term in our 



Lagrangian (2.1). In the broken phase we have an axion gradient in the near horizon 
geometry, which realizes a topological insulator interface, see Fig. [3] The character- 
istic patterning of the near horizon magnetic field, B = V x A, shown in Fig. |4| is 
reminiscent of the magnetoelectric effect at such interfaces. The magnetic vortices 
are localized near the black hole horizon and have alternating direction of magnetic 
field lines. 

In curved space the magnetic field is accompanied by vorticity, which is mani- 
fested by the function W. This causes frame dragging effects in the y direction. Test 
particles will be pushed along y with speeds W{r,x), in particular the direction of 
the flow reverses every half the period along x. The drag vanishes at the horizon 
and at the location of the nodes of W where x = n L/2, for integer n, see Fig. M In 
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Figure 2: Metric functions for T/T(. ~ 0.11. Note the metric functions A,B and C 
have half the period of W . The variation is maximal near the horizon, located at 
p = 0, and it decays as the conformal boundary is approached, when p — )■ cxd. 

general, the dragging effect remains bounded, and no ergoregion forms, where the 
vector dt becomes spacelike. 

3.2 The geometry 

There are several ways to envisage the geometry of our solutions, we discuss them in 
turn. 

The Ricci scalar of the RN solution is Rrn = —24, constant in r and independent 
of the parameters of the black hole. This is no longer true for the inhomogeneous 
phases, where the Ricci scalar becomes position dependent. Fig. [5] illustrates the 
spatial variation of the Ricci scalar, relative to the Rrn for T /Tf. ~ 0.054. The 
maximal curvature is always along the horizon at x = n L/2 for integer n. It grows 
when the temperature decreases and approaches the finite value of i? ~ —94 in the 



14 



CC 4- 

< 




-0.02 





2 0-2-4 0.5 1 



a. 




Figure 3: At relative to the corresponding RN solutions, Ay and ip for T/T^ ~ 0.11. 
The period of At is twice that of ip and Ay. The x-dependence dies off gradually as 
the conformal boundary is approached, at p — ;■ oo. 



small temperature limit. 

Embedding in a given background space is a convenient way to illustrate curved 
geometry. We consider the embedding of 2-dimensional spatial slices of constant x 



of the full geometry (2.2) 



as a surface in 3-dimensional AdS space 



dsl = S^!^dr^ + 2r2 e2^(^'-)dy2 



dsl 



dv 
2r^ 



(3.1) 



(3.2) 



We are looking for a hypersurface parametrized hy z = z{r). Then the metric on 
such a hypersurface reads 



dSn 



l + 2r^ 



dz 
dr 



dr 






c;2 j„,2 



2f 



- + 2r'dy 



(3.3) 
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Figure 4: Magnetic field lines for solution with T jTc — 0.07. The pattern of vortices 
of alternating field directions form at the horizon (located at p = 0). 



X 




Figure 5: Ricci scalar relative to that of RN black hole, R/Rrn — 1, Rrn = —24, 
for T/Tc ~ 0.054 over half the period. The scalar curvature is maximal along the 
horizon at x = n L/2 for integer n. 



Comparing (3.3) and (3.1) we obtain set of the relations 



2f 



l, + 2r-^^V 



dr J 



dr 



re 



c 



2r2/(0' 



(3.4) 
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Figure 6: The embedding diagram of constant x spatial slices, as a function of x at 
given y for T/Tc ~ 0.035. The geometry of p = const shces is maximally curved at 
X = n L/2 for integer n. 



resulting in the embedding equation 

^^Jf{r)-^e^B{r,.)-2Cir,x) _ (i + rdrC{r,x)y 



dz 

dr 2r^ 



(3.5) 



We integrate this equation for a given x, and in Fig. |6] show the embedding at 
constant y. The maximal curvature along p = const slices occurs at x = nL/2 for 
integer n, which is consistent with Fig. |5| 

The proper length of the stripe along x relative to the background AdS spacetime 
at given r is 

r.L/4 



lx{r)/lx{r = oo) 



(r, x) dx. 



(3.6) 



Fig. [7] shows the dependence of the normalized proper length on the radial distance 
from the horizon. The proper length tends to the coordinate length as 1/r^ asymp- 
totically as r — 7- oo, but it exceeds that as the horizon is approached. Namely, the 
inhomogeneous black brane "pushes space" around it along x, in a manner resembling 
the "Archimedes effect" . 



The proper length of the horizon in x direction is obtained calculating (3.6) at 



rg. Fig. M demonstrates the dependence of this quantity on the temperature. For 
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Figure 7: Radial dependence of the normalized proper length along x for T/Tc ~ 
0.054. While asymptotically the proper length coincides with the coordinate size 
of the strip, it grows as the horizon is approached. This is a manifestation of the 
"Archimedes effect" . 



high temperatures the length of the horizon resembles that of the homogeneous RN 
solution, however, it grows when temperature decreases. We find that at small T/T^ 
the proper length of the horizon diverges approximately as {T/Tc)~^'^. 

The transverse extent of the horizon, per unit coordinate length y, is given by 



ry[x) 



V2rne^'~'°'^\ 



(3.7) 



Fig. ^ shows the variation of ry{x) along the horizon for T/Tc ~ 0.054. Typically 
there is a "bulge" occurring at x = nL/2 and a "neck" at x = (2n + 1) L/4, for 
integer n. Comparing this with Fig. [5] we note that Ricci scalar curvature is maximal 
at the bulge and not at the neck as would happen, for instance, in the cylindrical 



geometry in black string case [31]. Fig. 10 displays the dependence of the sizes of 



the neck and bulge on T/Tc. Both sizes monotonically decrease with temperature, 
however the rate at which the neck is shrinking exceeds that of the bulge. This is 

In fact, we find that for ci = 4.5, r'^^^k j ^buige _ (T/Tcf/-i 
le range of temperatures that we investigated. For other 



demonstrated in Fig. 
near the lower end o 
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values of the axion coupling the scaling of the ratio is again power-law, with an 
exponent of the same order of magnitude, e.g. for ci = 8, the exponent is about 
0.12. This signals a pinch-off of the horizon in the limit T — )■ 0. 
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Figure 8: Temperature dependence of the proper length of the horizon along the 
stripe. Starting from as low as L at high temperatures, the proper length grows 
monotonically and for small T/Tc the growth is well approximated by the power-law 
dependence ~ (T/Tc 
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4 Thermodynamics at finite length 

In this section we consider the thermodynamics and phase transitions in the system, 
assuming that the stripe length is kept fixed. For the finite system the length of 
the interval is part of the specification of the ensemble and is kept fixed. In the 
next section we discuss the infinite system, for which the stripe width can adjust 
dynamically. 

4.1 The first law 

We demonstrated that below the critical temperature there exists a new branch of 
solutions which are spatially inhomogeneous. In the microcanonical ensemble the 
control variables of the field theory are the entropy S, the charge density N, and 
the length of the x-direction L, with corresponding conjugate variables temperature 
T, chemical potential fi, and tension in the x-direction tJq The usual first law is 
augmented by a term corresponding to expansions and contractions in the x-direction 
and is given by 

dM = TdS + lidN + tJL. (4.1) 

^Explicit expressions for these quantities in terms of our ansatz are given in appendix A 
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Figure 9: The the extent of the horizon in the transverse direction, Vy, as a functions 
of X for T/Tc ~ 0.054 in x G [—L/2,L/2]. The characteristic pattern of alternating 
"necks" and "bulges" forms along x. 
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Figure 10: The dependence of the size of the neck and the bulge on temperature. 
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where M, S, and N are quantities per unit length in the trivial y direction, but are 
integrated over the stripe. 



Our system has a scaling symmetry given by (2.19). In the field theory, this cor- 



responds to a change of energy scale. Under this transformation, the thermodynamic 
quantities scale as 



M^A^M, T^XT, S^XS, /x -> A/i, 



N ^ XN, 



Tx -^ A Tx, 



X 



Using these in (4.1) with A = 1 + e, for e small, yields 

2M = TS + fiN- T^L, 



(4.2) 



(4.3) 



the Smarr's-like expression that our solutions must satisfy and that can be used as 
a check of our numerics. For all of our solutions, we have verified that this identity 
is satisfied to one percent. 

4.2 Phase transitions 

The question of which solution dominates the thermodynamics depends on the en- 
semble considered. In the holographic context the choice of thermodynamic ensemble 
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is expressed through the choice of boundary conditions. The corresponding thermo- 
dynamic potential is computed as the on-shell bulk action, appropriately renormal- 
ized and with boundary terms rendering the variational problem well-defined. We 
examine each ensemble in turn. 

4.2.1 The grand canonical ensemble 

In our numerical approach, the natural ensemble to consider is the grand canonical 
ensemble, fixing the temperature T, the chemical potential fi, and the periodicity of 
the asymptotic x direction as L. The corresponding thermodynamic potential is the 
grand free energy density 

n{T,ij,L) = M-TS-fiN. (4.4) 

Different solutions of the bulk equations with the same values of T, fi, L correspond 
to different saddle point contributions to the partition function. The solution with 
smallest grand free energy Q is the dominant configuration, determining the thermo- 
dynamics in the fixed T, /z, L ensemble. In our case we have two solutions for each 
choice of T, /x, L, one homogeneous and one striped. Exactly how one one saddle 
point comes to dominate over the other at temperatures below the critical tempera- 
ture determines the order of the phase transition. 

In this ensemble it is convenient to measure all quantities in units of the fixed 
chemical potential fi. Then, after fixing L from the critical mode appearing at the 
highest Tc (see Fig. 18 and Table ^ in appendix B.l), we have that fi//!^ is a func- 



tion only of the dimensionless temperature T/jd. In the fixed chemical potential 
ensemble for large enough axion coupling we find a second order transition, where 
the inhomogeneous charge distribution starts dominating the thermodynamics im- 
mediately below the temperature at which the inhomogeneous instability develops. 
Near the critical temperature, the behavior of the grand free energy difference is 
consistent with {Q — Qrn)/^!'^ oc (1 — T/Tc)"^, while the entropy difference goes as 
{S — Srn)/^'' oc T/Tc — 1. This is as expected from a second order transition. As 



can be seen in Fig. [12] and Fig. [13} we find this second order transition for a range 
of lengths, L, and for a variety of values of the axion coupling Ci. With the current 
accuracy of our numerical procedure, we find it increasingly difficult to resolve the 
order of the phase transition for smaller values of ci. In fact, for ci = 4.5 the grand 
free energies of the homogeneous and inhomogeneous phases are nearly degenerate 
but still allow us to determine the phase transition as second order. It would be 
interesting to see if the phase transition remains of second order or changes to the 
first order for smaller values of the axion coupling. 
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Figure 12: The grand free energy relative to the RN solution for several solutions of 
different fixed lengths at ci = 8. In all cases shown we observe a second order phase 
transition. The critical exponents determined near the critical points in each case 
are consistent with the quadratic behavior {Vt — VLrn)/ l^^ oc (1 — T jTc)'^. 
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Figure 13: The grand free energy relative to the RN solution for c\ = 4.5 and fixed 
L/i/4 = 2.08. The grand free energies of the homogeneous and inhomogeneous phases 
are nearly degenerate, such that their maximal fractional difference is about 1%. 

To examine the observables in the striped phase further, we focus on ci = 8 and 
the corresponding dominant critical mode, Lfi/4 ~ 1.21, and consider solutions for 
the temperatures 0.00016 < T/Tc < 0.96. Various quantities are plotted with the 
corresponding homogeneous results in Fig. 14 Along this branch of solutions, the 



mass of the stripes is more than the RN solution and the entropy is always less. 
We plot the maximum of the boundary current density (jy), momentum density 
(Tj/o) and pseudoscalar operator vev (O^). Fitting the data near the critical point 
to the function (1 — T/Tc)"", we find the approximate critical exponents aj^ = 0.40, 
C(TyO = 0.41 and ao^ = 0.38 with relative fitting error of about 10%. 

We find evidence that the entropy of the striped black branes does not tend to 
zero in the small temperature limit, see Fig. [14| This is further supported by the 



behavior of the transverse size of the horizon (3.7). Here the bulge seems to contract 



at a much slower rate than the neck, which evidently shrinks to zero size in the limit 
T — !• 0. However, strictly speaking, this conclusion is based on extrapolation of the 
finite temperature data to T = 0. Checking whether the entropy asymptotes to a 
finite value or goes to zero in this limit, as suggested in [38], will require further 
investigation with a method of higher numerical accuracy. 

4.2.2 The canonical ensemble 

To study the system in the canonical ensemble we fix the temperature, total charge 
and length of the system. This describes the physical situation in which the system 
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Figure 14: The observables in the grand canonical ensemble for ci = 8 and I//i/4 = 
1.21 (points with dotted hne) plotted with the corresponding quantities for the RN 
black hole (solid line). Fitting the data near the critical point to the function (1 — 
T /Tc)^, we find the approximate critical exponents aj^ = 0.40, axyO = 0.41 and 
ao^ = 0.38 with relative fitting error of about 10%. 
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Figure 15: The difference in canonical free energy, at Ci = 8 and fixed length LN/A = 
1.25, between the striped solution and the RN black hole. The striped solution 
dominates immediately below the critical temperature, signalling a second order 
phase transition. 

is immersed in a heat bath consisting of uncharged particles. The relevant thermo- 
dynamic potential in this ensemble is the free energy density 



F{T, N,L) = M- TS. 



(4.5) 



If we measure all quantities in units of the fixed charge N, then, again, the free 
energy F/N"^ is only a function of the dimensionless temperature T/N. 

To solve our system with a fixed charge, we would need to fix the integral in x 
of the coefficient of the 1/r term in the asymptotic expansion of the gauge field Af. 
Numerically, it is much easier to fix the chemical potential, as this gives a Dirichlet 
condition on At at the boundary. In the grand canonical ensemble, we solved for one- 
parameter families of solutions at fixed L/i, labelled by the dimensionless temperature 
T/fi. Equivalently, in the {Lfi,T/fi) plane, we solve along the line of fixed L/x. 
Translated to the situation in which we measure quantities in terms of the charge 
density N, these solutions become one-parameter families of solutions with varying 
LN, or a curve in the {LN, T/N) plane with LN a function of T/N. By varying 
the length Lfi (or solving with /i = 1 and varying L), we can find a collection of 
solutions that intersect the desired fixed LN line. By interpolating these solutions 
and evaluating the interpolants at fixed LN, we can study the stripes in the canonical 
ensemble. 

In this ensemble we find a similar second order transition, in which the inho- 
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Figure 16: The entropy of the inhoniogeneous solution for Ci = 8 (points with dotted 
hne) and of the RN solution (solid line). Below the critical temperature, the striped 
solution has higher entropy than the RN. The RN branch terminates at the extremal 
RN black hole, while the striped solution persists to smaller energies. 



mogeneous solution dominates the thermodynamics below the critical temperature 



(Fig. 15). The scaling of the relative free energy density slightly below the critical 
temperature appears to fit a linear scaling, however using more points at lower tem- 
peratures in the fit increases the critical exponent towards a quadratic scaling 
expected in a second order transition. 



as 



4.2.3 The microcanonical ensemble 

The microcanonical ensemble describes an isolated system in which all conserved 
charges (in this case the mass and the charge) are fixed. This ensemble describes 
the physical situation relevant to the study of the real time dynamics of an isolated 
black brane at fixed length. In this case, the state that maximizes entropy is the 
dominant solution. As shown in Fig. [161 we find that the entropy of our inhonioge- 



neous solutions is always greater than that of the RN black hole of the same mass. 
Furthermore, the mass of the inhomogeneous solutions is always smaller than that 
of the critical RN black hole. Therefore, at fixed LN , the unstable RN black holes 
below critical temperature are expected to decay smoothly to our inhomogeneous 
solution. 
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4.2.4 Fixing the tension 

Alternatively, one could attempt to compare solutions with different values of L. 
The meaningful comparison is in an ensemble fixing the tension Tx- For example, 
one could compare the Legendre transformed grand free energy 

G{T,^^,Tx)=M-TS-^lN-TxL (4.6) 

where the additional terms comes from boundary terms in the action rendering the 
new variational problem (fixing Tx) well-defined. The candidate saddle points are the 
solutions we find with various periodicities L, and their relative importance in the 
thermodynamic limit is determined by G{T,ij,,Tx)- In particular the solution which 
is thermodynamically dominant depends on the value of Tx we hold fixed. In this 
study we concentrate on the thermodynamics in the fixed L ensemble and we leave 
the study of the fixed r^ ensemble to future work. 

5 Thermodynamics for the infinite system 

In this section we lift the assumption of the finite extent of the system in the x- 
direction and consider the thermodynamics of the formation of the stripes below the 
critical temperature. For the infinite system we can define densities of thermody- 
namic quantities along x: 

MSN 
m=— , s = -, n=j^- (5-1) 

In terms of these, the first law for the system becomes 

dm = Tds + fj,dn (5.2) 

and the conformal identity is 

3m = 2{Ts + fin). (5.3) 

In the infinite system, we compare stripes of different lengths, at fixed T/fi, 
to each other and to the homogeneous solution. The solution that dominates the 
thermodynamics is the one with the smallest free energy density u, where 

u = m — Ts ~ fin. (5-4) 



This comparison is shown in Fig. 17 for ci = 8, where we see that the free energy 



density of the stripes is negative relative to the RN black hole, indicating that the 
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Figure 17: Action density for ci = 8 system relative to the RN solution. The red 
line denotes the approximate line of minimum free energy. 



striped phase is preferred at every temperature below the critical temperaturejj Very 
close to the critical temperature, the dominant stripe is that with the critical wave- 
length kc- As we lower the temperature, the minimum of the free energy density 
traces out a curve in the (L/x, T) plane, and the dominant stripe width increases to 
L/i/4 ^ 2. 

One can also study the observables of the system along this line of minimum free 
energy density. The results are qualitatively similar to those for the fixed L system 



(Fig. 14). In particular, the free energy density scales as {uj — urn)/!^ cx: (1 — T/Tc 



near the critical point, indicating a second order transition in the infinite system as 
well. 
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^In appendix B.4 we describe the generation of Fig. 17 
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A Asymptotic Charges 

A.l Deriving the charges 

Since our ansatz is inhomogeneous and includes off-diagonal terms in the metric, 
and our action is not standard (in that it includes the axion coupling) we have re- 
derived the expressions for the charges and other observables in our geometry. In 
deriving the asymptotic charges of our spacetime, for the four dimensional Einstein- 
Maxwell- Higgs theory we discuss in the main text, we follow the covariant treatment 
of |ini HI] • We refer the reader to those papers for details of the method used. 

The bulk action has to be supplemented by boundary terms of two types. First, 
there are boundary terms needed to ensure that the variational problem is well- 
defined. Then there are counter-terms, terms depending only on the boundary values 
(leading non-normalizable modes) of fields on the cutoff surface, which are added to 
render the on-shell action and the conserved charges finite. Both kinds of boundary 
terms are the standard ones for Einstein-Maxwell-Higgs theory; the additional axion 
coupling does not necessitate an additional boundary terms of either kind as long as 
the scalar mass satisfies m^ < 0. 

We find it convenient to study the first variation of the on-shell action, which 
always reduces to boundary terms. The expression for the regulated first variation 
of the on-shell action can be differentiated with respect to the boundary values of 
the bulk fields, to give finite expressions for the conserved charges. We write those 
expressions below in terms of the asymptotic expansion of the fields occurring in our 
ansatz, carefully taking into account the differences between our coordinate system 
and the standard Fefferman- Graham form of the asymptotic metric, which is used 
to derive the standard expressions in the literature. 

Having explained our procedure, we now display the expressions for the observ- 
ables used in the main text. We first assume the radial coordinate is in the standard 
Fefferman- Graham form, and then discuss additional terms arising from change of 
coordinate necessary to bring our asymptotic metric into the standard form. 

For the scalar fields ip, one can write asymptotically 

ij{x,r) =^W(x)r-^- +ij^^\x)r-^+ (A.l) 
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with 

X± = l± /f+^- (A.2) 

We set ip^'^^x) = as part of our boundary conditions, then the coefficient ip^^\x) 
is the spatially modulated VEV of the scalar operator dual to ip. 
Similarly, the gauge field can be expanded near the boundary as 

A,{x,r) = A'^^\x)-^^. (A.3) 

The functions Afj, \x) correspond to the charge and current density of the boundary 
theory. 

As for the boundary energy-momentum tensor, the expression is fairly simple in 
odd number of boundary dimensions, and we have checked that it is not modified by 
the matter action. With our normalization convention one can write 

T., = Gglf, (A.4) 

where the superscripts of the metric functions denote the order in the asymptotic 
expansion. 

Since our metric ansatz is not of the Fefferman- Graham form, we need to perform 
a change of coordinate (in the x, r plane, for which we used the conformal ansatz) to 
put the metric is such a form. The details of the transformation are straightforward 
and the process results in the following shifts in the asymptotic metric quantities: 



,(3) 



2 



A^^^ = -9{x), (A.5) 



for every i, j, where g{x) is the leading asymptotic correction to the metric component 
Qrr- That is, at large r that metric component becomes 



]_ g{x) 
2r^ r^ 



9rrir,x)^— + ^. (A.6) 



Finally, since the metric becomes diagonal asymptotically, the non- vanishing time 
components of the energy-momentum tensor T^ and Tyt have a simple interpretation 
as energy and momentum density, respectively. The conserved charges are given by 
integrating those densities over a spatial slice. 
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A. 2 Explicit expressions for the charges 

A. 2.1 Homogeneous solution 

For reference, in tliis subsection we give the explicit expressions for the homogeneous 
RN solution in our conventions. The radius of the horizon is given in terms of the 
temperature by 

ro = ^ (2nT + ^3fi^ + An^T^^ . (A. 7) 

The mass, entropy and charge of the RN solution of fixed length L are 



Mrn = (4ro' + Ao)^, (A.8) 

Srn = ^TrrlL, (A.9) 

Nrn = 2ro/iL. (A. 10) 

The corresponding densities in the infinite system are given by dividing through by 
L. 

A. 2. 2 Inhomogeneous solution 

Here we list explicit expressions for the thermodynamic quantities in our system in 
terms of our solution ansatz. Conserved charges are given by integrating over the 
inhomogeneous direction . W e define /'•^•' = — (Irj] + /i^ro)/4, the 1/r^ term from the 



function /(r) (equation (2.3)), and X^^\x), for X = {R,S,T}, as the coefficient of 
the 1/r^ term of the corresponding metric function. The energy-momentum tensor 
yields the mass 
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M= [ {T'\i))dx = 4 /" ^{xf{-f^^^ + hS^^\x) + ?,T^^\x))dx, (A.ll) 

Jo Jo 

the tension in the x direction 

T^ = - f (T^^(x))rf5 = 2 f ^{xy{f^^^ + 6R^^\x)+AS^^\x) + 6T^^\x))dx, (A.12) 
Jo Jo 

and the pressure in the y direction 

Py= {Tyy{£))dx = -2 / ^(a;)=^(/(^) + QR^^\x) + 10S^^\x))dx. (A.13) 

Jo Jo 



^''See appendix B.3 for details about the numerical process, including the definitions of the x 
coordinate and £,{x). The functions {R,S,T} are defined on the UV grid; they are analogous to 
{A, B, C} in the original ansatz. 
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Now, expanding the equations of motion at the asymptotic boundary, we get the 
relation R^^\x) + 2S^^\x) +T^^\x) = 0. Using this, we see that {T'"'{z)) is traceless, 
as necessary. Conservation of the energy momentum tensor requires dxTx = 0. This 



is related to the constraint equation (2.18) and we explain our strategy to ensure it 



is satisfied in appendix B.3 



The coefficient of the 1/r falloff of the gauge field gives the charge 

N = -2 aI^\x). (A.14) 



At the horizon, we read the (constant) temperature as 

T = —il2rl - /i2)e-(^-^)l'-='-o (A.15) 

and the entropy is proportional to the area of the event horizon, given by 

S = 47Trl / e(^(''«'")+^(''0'"))da;. (A.16) 

Jo 

A. 3 Consistency of the first laws 

Here, we discuss the first laws for both the finite length stripe and the stripe on the 
infinite domain. 

A. 3.1 Finite system 



In our system, as described above, we have unequal bulk stresses r^^^ and Py. Then, 
if we have a rectangle of side lengths {L,Ly), the work done by the expansion or 
compression of this region will differ depending on which direction the stress is in. 
The usual —PdV term in the first law is replaced and we have 

dM = TdS + fidN + TxLydL - PyLdLy, (A.17) 

where the hatted variables represent thermodynamic quantities integrated over the 
entire system. Defining densities (in the trivial y-direction) by 

M^f S^l N^f, (A.18) 

1^1, J-^ll -^11 



^We define r^ = —Px^ where Px is the pressure in the x direction. For our solutions, Tx > 0. 
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we can write the first law as 



dM = TdS + fidN + T^dL + —^{-M + TS + fxN- PyL). (A.19) 

Ly 

Tracelessness of tlie energy-momentum tensor implies M = L{Px + Py), so that the 
term proportional to dLy can be rewritten as the conformal identity (4.3), which 
disappears for a conformal system described by the first law (4.1). Therefore, the 



first law (4.1) and the conformal identity (4.3) are consistent 



A. 3. 2 Infinite system 

For the infinite system, we define densities in both the x and y directions as equation 



(5.1). Under the scaling symmetry (2.19), these scale as 

m — )■ A^TTi, s — )■ A^s, n — )■ A^n. 



(A.20) 



Using the first law (5.2), we derive the conformal identity (5.3). Again, we can see 



this from the first law for the system with integrated charges. Plugging the densities 



m, s,n into the first law of the finite length system (4.1), we arrive at 



dm = Tds + fidn -\ — — (— m + Ts + ^s + r^) 



(A.21) 



Using the conformal identity of the finite length system (4.3), we see that the term 



proportional to dL is just the conformal identity for the infinite system, which is 



satisfied for a system described by (5.2). 



B Further details about the numerics 

B.l The linearized analysis 

Following [22], we look for static normalizable modes around the Reissner- Nordstrom 
background. We consider the fiuctuatioE 



12 



Sgty 
6A. 



A 



(r - ro) 



w{r) sin{kx] 



A(a(r) sin(/cx)), 
A(0(r) cos{kx)), 

^^Regularity at the black hole horizon enforces that Sgty{ro) — 0. 



y 

6ij 



fB.l) 
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where A is a small parameter in which we can expand the equations. Putting this 



ansatz into (B.6) - (B.12) and expanding to linear order in A, we arrive at the 



linearized system 

/// fo(i'(r) (4r — 2ro)w' (r) w(r) (2rQ (4r^ + At'^tq + Arro^ — tq) — k'^r'^) 

w (r) 1 1 : = 

r^(r — ro) r(r — tq) r^ (4r'^ — r (4ro'^ + tq) + tq^) 

^ (8r^ + r {Arp^ + rp) - 2ro^) a'{r) _ 

r (4r4 - r {Atq^ + Tq) + Tq^) 
k'^a{r) Cikro(f){r) 



4rro'u;'(r) 4ro^w(r) 



4r^ _l_ 47-2^Q _|_ 47-7-q2 _ j,^ 4^4 _ ^ (4ro^ + tq) + ro^ 

2y3r2 (4r4 - r (4ro3 + tq) + tq^) 
0(r) {k^ + 2mV2) (-16r^ + 4ro^ + tq) 0'(r) 

4^4 _ J. (47-q3 _|_ ^^^ _|_ ^^2 4^4 _ j^ (4ro^ + To) + ro^ 



0, (B.2) 



0. 



Fixing the scalar field mass as m'^ = —4, there are three parameters in these equa- 
tions: the temperature of the black brane To (equivalently the location of the horizon 
To), the wavenumber k, and the strength of the axion coupling ci. In this analysis, 
we will choose ci and k and then use a shooting method to find the Tq at which 
normalizable modes appear. 

Due to the linearity of the equations, the scale of our solutions is arbitrary. We 
use this to fix a Dirichlet condition on w at the horizon. Changing coordinates to 
P = v^r^ — Tq, and expanding the equations near p = gives regularity conditions 
on the fiuctuations at the horizon in terms of Neumann boundary conditions. Our 
horizon boundary conditions are then 

wip)l=o = 1, w'ip)l=o = a'{p)l=o = 0'(p)|p=o = 0, (B.3) 

Namely, that the fields are quadratic in p near the horizon. In order to search for 
normalizable modes, we set the sources in the field theory to zero by imposing leading 
order fall-off conditions near the AdS boundary: 

^(p) = ^ + ..., a(p) = ^ + ..., 0(p) = ^ + .... (B.4) 

p-^ p p^ 
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Figure 18: The critical temperatures at which the Reissner Nordstrom black brane 
becomes unstable, for varying axion coupling ci. As the strength of the axion cou- 
pling increases, the size of the unstable region (the area under the critical temperature 
curve) also increases. 

In practice, after fixing ci and k, we use Tq as a shooting parameter to find the 
solution with the correct w fall-off and the corresponding critical temperature T^.. 

For each Ci, we find a range of unstable momenta. By adjusting the strength of 
the axion coupling, one can find a large variation in the size of this unstable region in 



the (A;//i, To//i) plane (see Fig. 18). The relationship between ci and the maximum 
critical temperature is well fit by T™"^(ci)//i = 0.025ci — 0.091. The wavenumbers 
for the dominant critical modes, corresponding to r™"'*(ci), for select Ci are found 
in Table [U 



Cl 


rp^axl^ 


kJii 


L/i/4 = 7r/2fcc 


4.5 


0.012 


0.75 


2.08 


5.5 


0.037 


0.92 


1.71 


8 


0.11 


1.3 


1.21 


18 


0.37 


2.85 


0.55 


36 


0.80 


5.65 


0.28 



Table 1: The maximum critical temperatures and corresponding critical wavenum- 
bers for varying Ci. 
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B.2 The equations of motion 

For completeness, here we present the equations of motion derived from the La- 



grangian (2.1). The Einstein equations in our case are four second order eUiptic 
equations for the metric components and two constraint equations. For the com- 
pactness of the expressions, we define 

du-dv = drUdrV + ^d,ud,v, o^u = dlu + -^.dlu. (b.s) 

The four elhptic equations, formed from combinations of G\ — T^ = 0, Gy — Ty = 0, 
Qy - Ty = 0, and G; + G^ - (T; + Tf ) = 0, then take the form 

6^ A + {OAf + dA-dc- "" {OWf - "^^AOAtf 



and 



1 fe-^^W^ 2c^./^.^2 e-2^W^^, ^, /5 3/'\ , 



.1 f'\^^ 3 3e2^ e^^mV g^/ f. 



-1 -2A+2C -| 

O^B + -(Otpf - , , (OWf -OA-OC- -drA 
2 4/ r 



— 2AH-2C* —1A 

d^G + (6c)2 + OA ■ 6c + "^ (oiy)^ + ^(OA)^ 

+ l5A+/^^+n5C'+-- — + ^'''^'^' + ^ = (B 8) 

fr- \ T f / V T f A/r^ f f f 



O'W -OA-OW + WG ■ OW - ^ ^ {OAyf 

e-2C . . 4 
—OAfOAy + -drW = Q. (B.9) 
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The matter field equations are 

—A—C 

O^tp + dA-d^ + dC-d^+ '^^ {drAtd^Ay - d^AtdrAy) 

8\/3r^f 



O^At - OA ■ OAt + OC ■ OAt + OW ■ OAt + OW ■ OA^ 

f 

-2A+2Ctx^2 

+2W0C ■ OAy - 2W0A ■ OAy + OW ■ OAy 



4x/3r2 V / 



+ ^^ ( e-^-^ ) {dri^d^Ay - d^i^drAy) 



idr^d^At-d^lPdrAt) + -drAt-^drAy = Q, (B.ll) 



Ay/Sr^f r f 



and 



~2A+2C^r -2A+2C 

6^ Ay + 6 A ■ OAy - OC ■ OAy " OW ■ OAy — OW ■ OAt 



+ ! /7^ o r (dri'd.^At - d,^^dr.At) + \ /- ^ , {dr^pd^Ay - d^i)drAy) 
4v3r^/ Ay6r^j 

+ f - + 4 ) drAy = 0. (B.12) 



2 /'~ 

- + — 
r f 

Finally, the constraint equations are 

-2e~^^ fr\d^drA + d^drC) + 2e~^^ fr^drA {d^B - d^A) 

+2e-^^fr^ {d,A + d,C) drB + 2e-^^ fr^ {d^B - d,C) drC - e-^^r^f'd,A 

^^2g-2(A+i?-c)^^^^^^ _ /e-2(^+^)9,A,9,A, - 2e-2^/r29,^9,^ = 

(B.13) 
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and 



^l^ + die - ^MA + dlO + (l - ^) {drAf + (l - 4^) {drCf 



f f 2\ P"2A -2A 



2 



^-''^^-^^ A. „,,2 1 ,. ...^2\ , e-2^ A, , ,2 1 .o . ^2 



B.3 Constraints 

The constraint equations, Gl^—T^ = and C^—G^ — iT'^—T^) = 0, are the non-trivial 
Einstein equations that are not part of the system of second-order elhptic equations 
that we numerically solve. As discussed in ^ the weighted constraints can be shown 
to solve Laplace equations on the domain. If we satisfy one of the constraints on 
all boundaries and the other at one point, they will be satisfied everywhere. At the 
black hole horizon, we choose to impose t-^\/7\/— ^(GJ! — G^ — (T^ — T^)) = at the 
point (p, x) = (0, 0) and a/^(G^ — T^) = across the horizon. Since we use periodic 
boundary conditions in the inhomogeneous direction, the boundaries at x = and 
X = Xmax are trivial if y/—g{Gl. — T^) = at the horizon and the conformal boundary. 
Then, we are left with the task of satisfying ^—giG"^^ — T^) = at the boundary. 
In ® we found the asymptotic expansion of this constraint as 



39,A(3)(x) + 29,S(3)(x) + 39,C(3)(x 



G: - T: (X -^ ' ' ' ,' ' ' ^ ^ + 0{r~% (B.15) 



.2 



where A^^\x) ^ B^^\x) and C^^\x) come from solving the elliptic equations. It ap- 
pears that, within our problem, we do not have the ability to make the weighted 
constraint disappear. The key lies in an unfixed gauge symmetry in our original 



metric that is related to conformal transformations of the (r, x) plane/*^^ Essen- 
tially, within our metric ansatz, we have the freedom to transform to any plane 
(r',a;') that is conformally related to (r, x). Demanding that the weighted constraint 
sJ—giG"^^ — T^) vanishes at the conformal boundary uniquely identifies the correct 
coordinates (f, x). 

^■^See [33] for a discussion of the same issue in a different context. 
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Our procedure is to split the domain at some intermediate radial value pint- On 
the IR portion of the grid, < p < pint-, the equations are as above. On the UV 
portion of the grid, pmt < p < Pcut, we use the coordinate freedom to select the 
correct asymptotic radial coordinate. We can write the metric in the UV as 

^^2 ^ _2r~2/>,x)e2«rft2 + e^^ ( ^ff ^ + 2f^dxA + 2f^e'^{dy - Udt)\{B.16) 

\Zr J (r, X) J 

where f{r,x) = f{r{f,x)). Under a transformation in the (f,x) plane such that f 
and X satisfy Cauchy-Riemann-like relations 

dr{r,x) f{r,xY dx{r,x) dx{r,x) 1 dr{r,x) 



Or r^ dx ' Or 4r'^r[r,xyf{r) dx 

the metric becomes 



(B.17) 



ds^ = -2Hr,x?f(r)e'^^dt^ + e'^^\Vr(r,x)\^ ( -^tttt + 2r^dx^ ] + 2Hr,xYe^^(dy - Udtf 

2r^j[r) 



{B.li 



|Vr>,x)P = ^rg%^V+ .,., \.... f^%^V- (B.19) 



with 



r{r\xY \ dr J Ar'^r(r,x)'^f{r) \ dx 

We now have an extra function f{r,x) in our system which we may use to satisfy 
the constraint and fix the residual gauge freedom, as we will now see. The Cauchy- 
Riemann-like conditions give the Laplace-like equation 

or \r(r, xj^ or j ox \Ar'^r[r,x)^j{r) ox J 
We can solve this asymptotically, finding 

r[r, x) = a^)r + ^ 24e(x)r + ' ' ' ' ^^'^^^ 

where C,{x) is an arbitrary function that encodes the coordinate freedom we have. 
Expanding the constraint asymptotically, we have 



Gl - T; oc \ (2{3d^R^^\x) + 2d^S^^\x) + 3d^T^^\x))^{x) 

+ 3(/(3) + 2R^^\x) - 4S^^\x) + 2T^'^\x))^'{x)') + 0{r-^), 

(B.22) 
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where X = X^^\x)/r^ + . . . asymptotically, for X = {R,S,T}. Demanding that 



the constraint (B.22) vanishes at the leading order yields a differential equation we 
can solve for C,{x), giving us a boundary condition for the function r{r,x), such that 
the weighted constraint will disappear at the conformal boundary. However, we have 
found that the code is unstable if we directly use this solution for ^(x). Instead of 
directly integrating the constraint, we use the freedom in ^{x) to fix the tension r^, 
to be constant. This enforces the same effect on the tension as if we had used the 
explicit solution for ^(x) but is much more stable numerically. Below, we check that 
the constraints are suitably satisfied even though our boundary conditions do not 
exactly fix them. To this end, we set 



ax) 



K 



(/(3) + 6i?(3)(x) + 45(3) (x) + 6T(3)(x))i/3- 



(B.23) 



Expanding the equations asymptotically gives the expression R'^'^\x) + 2S^'^\x) + 
T^^\x) = 0; if this is satisfied on our solutions our definition of ^{x) coincides with 



that found by integrating the constraint (B.22) 



The constant K appearing in ^{x) sets the scale of the boundary theory. We use 
it to fix the length of the inhomogeneous direction in the field theory to be L/i/4. 
The correct coordinate in the inhomogeneous direction of the field theory is x. From 
the Cauchy-Riemann conditions, we can find the large r expansion of x(r, x) as 



x[r,x] 



^ dx' 



+ 



f(x) 



a^') 8e(x)2r 



+ .... 



(B.24) 



Integrating to find the proper length of one cycle in the boundary, we solve for K at 
leading order in r to find 



K 



4 

z 



L/4 



(/(3) + m^''\x) + AS^^\x) + QT^^\x)f/^dx'. 



(B.25) 



When integrating the charges over the inhomogeneous direction in the field theory, 
one must remember to integrate over the correct coordinate, dx = dx/^{x). 

Our corrected numerical procedure is as follows. On the IR grid, we solve the 
elliptic equations (B.6) - (B.12) for the metric functions A,B,C and W. On the 



UV grid, we solve the equivalent elliptic equations from the metric (B.18) in the 



variables i?, S", T and U plus equation (B.20) for the new field f{r\ x). At the horizon, 
we enforce the boundary conditions discussed in ^ At the interface p = pmt, we 
impose matching conditions on the four metric functions and that r{pint, x) = r{pint)- 
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=1, 




Figure 19: The data underlying Fig. 17 The points represent solutions we computed. 
These were interpolated to find the free energy density over the domain. The solid 
blue line is the edge of the unstable region and the thick red line is the approximate 
line of minimum free energy density. 



Asymptotically, R,S,T and U all fall off as 1/r^. To set boundary conditions on r, 
we notice that 

^3 



drr{r,x) + 



r 



2e(x) + O 



r- 



(B.26) 



We truncate this expression at 0{r^ ) and finite difference to find an update proce- 
dure for r{pcut,x). This boundary condition is updated iteratively as the functions 
R, S, T are updated in our solving procedure such that once we find a solution with 
small residuals we can be sure that the tension is constant and the constraint is 
satisfied. 



B.4 Generating the action density plot 

To generate the relative action density plot, Fig. [17], we find the solutions on a grid of 
lengths L and temperatures Tq, as shown in Fig. 19 By interpolating these solutions 
on the domain, we can map the thermodynamic quantities across the unstable region 
and determine the approximate line of minimum free energy, or the dominant solution 
in the infinite size system. 
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B.5 Convergence and independence of numerical parame- 
ters 

B.5.1 Performance of the method and convergence of physical data 

As discussed above, to solve the equations numerically, we use a second order finite 
differencing approximation (FDA) before using a point-wise Gauss-Seidel relaxation 
method on the resulting algebraic equations. The method, including the UV proce- 
dure described above, performs well for this system. 

The UV procedure is unstable for a generic initial guess, resulting in a divergent 
norm. To find a solution from a generic initial guess, we can run the relaxation 
without the UV procedure until the norm is small enough that the result approxi- 
mates the true solution, before activating the UV procedure to find the true solution. 
Once we have these first solutions, by using these as an initial guess for solutions 
nearby in parameter space and by interpolating to a finer grid, we can generate fur- 



ther solutions by relaxing with the UV procedure. In Fig. 20, we plot the L^ norm 
of the total residual during the relaxation of the ci = 8 solution at Tq = 0.04 and 
L///4 = 0.75 for the grid spacings dp,dx = 0.04,0.02,0.01, showing the expected 
exponential behavior of the Gauss-Siedel relaxation. The physical data extracted 
from our solutions is consistent with the expected second order convergence of our 



FDA scheme, see Fig. 21 



B.5. 2 Asymptotic versus first law mass 

A useful check of the numerics is to compare the mass of the system read off from 



the asymptotics of the metric, equation (A. 11), to that computed by integrating the 



first law, equation (4.1). Since the temperature and entropy are read off from the 
horizon, comparing these two methods of finding the mass provides a non-trivial 
global consistency check on our results. We verify that the difference between the 
asymptotic mass and the first law mass remains smaller than 0.5% across our set of 
trials, indicating consistency of our results. 

A related check of the numerics is the conformal identity or the Smarr-like rela- 
tion, 2M = TS + fiN — T^L, derived above from the first law for the finite length 
system. To evaluate how well our solutions satisfy this equation, we examine the 

2Mfau^off-TS-fiN + r^L ^ 
meix{Mfaii-off,TS,nN,TxL)' 

since the largest term in the expression sets a scale for the cancellation we expect. 
This ratio is very small for our solutions near the critical temperature. As we lower 
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Figure 20: The behavior of the L^ norm of the residual during the relaxation itera- 
tions for Ci = 8, To = 0.04 and L/i/4 = 0.75. From top to bottom (at the left of the 
plot) the grid spacing is dp, dx = 0.04, 0.02, 0.01. The UV procedure is unstable un- 
less the solution is close enough to correct solution. For grid spacing dp, dx = 0.04, 
the UV procedure was activated after 3 x 10^ iterations while for the others, the 
initial guess was taken to be a solution with slightly different parameters such that 
the UV procedure could be used immediately. 

the temperature, this ratio increases, but stays small. The precise value depends on 
the parameters of the solution, but is not larger than order 1%. Moreover, this ratio 
decreases as we move the position of the finite cutoff of the conformal boundary to 
a larger radius. 

B.5.3 Finite pcut boundary check 

For the ci = 8 trials reported in the paper, we use pcut = 12 as our conformal 
boundary. In Table |2] we present results for varying p^ut, showing that our choice is 
large enough such that the physical results are insensitive to the cutoff. Although the 
physical results presented in the table appear very stable, at small pcut, the results for 
the mass and charge depend significantly on the fitting procedure for the asymptotic 
metric functions and gauge field. By running our simulations at pcut = 12, we are 
both well within the the region where the solutions do not change with the conformal 
boundary and within a region where our fitting procedure to the asymptotics behaves 
well. 
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Figure 21: The value of the scalar field condensate for varying grid sizes for ci = 8 
and Lfj,/4: = 0.75. From top to bottom, the grid spacing is dp, dx = 0.01, 0.02, 0.04. 
The results are consistent with second order scaling as expected from our numerical 
approach. 

B.5.4 Behavior of the constraints 

One of the most important checks for our numerical solution is the behavior of the 
constraints. For numerical homogeneous solutions found with our method, the L^ 
norm of the constraints is very small, on the order of 10~^. For the inhomogeneous 
solutions, the constraints are small near the critical temperature, but grow and sat- 
urate as we lower to the temperature, to have a maximum L^ norm on the order 



of 10~^: see Fig. 22 Since our boundary conditions explicitly fix the weighted con- 
straints on the horizon, they disappear there. The weighted constraints then increase 
towards the conformal boundary, approaching a modulated profile of constant am- 
plitude. The amplitude near the conformal boundary controls the overall L^ norm 
of the constraints. 

The constraint violation improves marginally with step size and with moving 
the interface closer to the horizon, but does not improve as we take the conformal 
boundary to a larger radius. To check that the constraints are well satisfied on our 
solution, we compare them to the sum of the absolute value of the terms that make 
up the constraints. That is, if the constraints are given by Yli ^ii '^^ compare this 
to Ylii \^i\- This procedure gives us an idea of the scale of the cancellation among 
the individual terms hi. We find that the sum ^^ \hi\ diverges approximately as r^ 
towards the asymptotic boundary, such that the approach of the constraint violation 
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^(G^ - Tr, 



Wv^(g;-g^-(t;-t;)) 



0.005 

0.000' 
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Figure 22: The weighted constraints for ci = 8 and Lfi/A = 1.21. The top plots are 
near the critical point, T/T^ = 0.97, while the bottom plots are at small temperature, 
T/Tc = 0.00016. By our boundary conditions, the constraints disappear at the 
horizon. They approach a finite value as they approach the asymptotic boundary. 
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Pcut 


S 


M 


A^ 


1 


0.758504 


0.305774 


0.527406 


2 


0.767913 


0.342327 


0.490524 


3 


0.768211 


0.341928 


0.490593 


4 


0.768285 


0.342043 


0.490583 


5 


0.768311 


0.342136 


0.490577 


6 


0.768322 


0.34221 


0.490574 


7 


0.768328 


0.342277 


0.490572 


8 


0.768332 


0.342324 


0.49057 


9 


0.768334 


0.342367 


0.490569 


10 


0.768335 


0.342402 


0.490568 


11 


0.768336 


0.342434 


0.490568 


12 


0.768336 


0.342459 


0.490567 



Table 2: Behavior of physical quantities with the cutoff for ci = 8 and Lfi/A = 0.75 
and for fixed grid resolution dp, dx ~ 0.02. The entropy S is read off at the horizon, 
while the mass M and the charge A^ are read off at the conformal boundary. Both 
the entropy and the charge are very robust against the location of the conformal 
boundary. The mass takes slightly longer to settle down, but is well within the 
convergent range for pcut = 12. 



to a constant is a good indicator that the constraints are satisfied on the solution. In 
Table [3| we compare the L^ norm of these two sums on the entire domain, showing 
that the constraint violation for the inhomogeneous solutions is generally about four 
orders of magnitude less than the scale set by ^^ \hi\. Interestingly, the relative 
constraint improves marginally as we go to lower temperatures. 



B.5.5 The asymptotic equation of motion 

Expanding the equations of motion asymptotically gives the relation 



(B.28) 



which can be used to give another check of the numerics. As explained in A. 2. 2 
this condition implies the tracelessness of the energy- momentum tensor. For the 
inhomogeneous solutions near the critical temperature we find that this expression 
is on the order of the individual metric functions X'^^\ where X = {R,S,T}, but 
generally decreases as we lower the temperature. As well, we find that homogeneous 



solutions found using our numerical techniques satisfy (B.28) well. There seems to 
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Parameters 


To 


L\Y..hi)lL\Y.. 


h^) 


ci = 8, L/i/4 = 


= 2.00 (RN solution) 


0.105 


9.12-10-^ 




Ci = 8, L/i/4 = 


= 1.21 (striped solution) 


0.075 
0.05 
0.025 
0.005 
0.001 


2.02 ■ 10"^ 
1.84 ■ 10"^ 
1.58 ■ 10-^ 
1.37-10-^ 
1.32 ■ 10-^ 





Table 3: Comparison of the constraint violation, measured by the schematic con- 
straint equation J2i^ij ^o ^^e scale set by the individual terms, XIjI^jI; ^^ g"^ 
size dp, dx ~ 0.01. We take the L^ norm of the measures on the entire domain. 
The Ci = 8, L/i/4 = 2.00 solution is a homogeneous RN solution found numerically 
with our code, for which the constraints are very well satisfied. The constraints for 
the striped solutions are satisfied compared to the scale set by ^^ \hi\ by four or- 
ders of magnitude and the relative constraint improves marginally as we lower the 
temperature. 

be an unidentified systematic error here that may deserve further attention in the 
future. Possible problems may occur in the implementation of the UV procedure or 
in our procedure to read off the coefficients of the falloffs of the metric functions. 
However, our physical results are robust under changes to the boundary conditions, 
so that we are confident in our results despite this possible systematic. In particular, 
the physical quantities extracted from the horizon are independent of the different 
boundary constraint fixing schemes we implemented. Therefore, we advocate using 
the mass derived from the integrated first law, which uses no asymptotic metric 
functions. 
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